import numpy as np
import matplotlib.pyplot as plt
import pandas as pd

import igrins_mod as ig # Custom module file for igrins shenanigans

import glob
import os
import warnings

import astropy.units as u
from astroquery.nist import Nist # atomic lines
# from astroquery.linelists.cdms import CDMS # molecular lines?

from lmfit import Model, Parameters, CompositeModel
from lmfit.models import GaussianModel

from scipy.stats import chisquare
from scipy.integrate import trapz, simpson
from scipy.optimize import curve_fit

from astropy.io import fits

# Plotting Parameters
plt.rcParams['font.size'] = 20
plt.rcParams['axes.labelsize'] = 18
plt.rcParams['xtick.labelsize'] = 18
plt.rcParams['ytick.labelsize'] =18

plt.rcParams['legend.fontsize'] = 16
plt.rcParams['figure.titlesize'] = 20

plt.rcParams['axes.labelweight']='bold'
plt.rcParams['axes.linewidth'] = 3

plt.rcParams['xtick.major.size'] = 10
plt.rcParams['xtick.minor.visible'] = True
plt.rcParams['xtick.minor.size'] = 5

plt.rcParams['ytick.major.size'] = 10
plt.rcParams['ytick.minor.visible'] = True
plt.rcParams['ytick.minor.size'] = 5

plt.rcParams['ytick.direction'] = 'in'
plt.rcParams['xtick.direction'] = 'in'

%matplotlib inline

Things that should only have to be defined once, here at the top

---

# Size of 1 spectral resolution element
# IGRINS Spectral Resolution
spec_res = 0.00001

c = 299792458  # speed of light m/s

# Use Normalized (single) Gaussian Distribution
def gaussian_func(x, ampl, center, std):
    return ((ampl) / (std * np.sqrt(2 * np.pi)) * np.exp(-0.5 * ((x - center) / std) ** 2)) + 1

# Reduced and order-merged data filepath 
# Laptop Path
data_path = "C:\\Users\\Savio\\Documents\\GitHub\\IGRINS-Spectra\\IGRINS_Merged"

# File path for figures to live in
fig_path = "C:\\Users\\Savio\\Documents\\GitHub\\IGRINS-SpectraIGRINS_figs\\standards_spectra"

# Create the folder if it doesn't exist
if not os.path.exists(fig_path):
    os.makedirs(fig_path)

# Nicole's merged K-band spectra of some Taurus Standards
merged_standard_files = glob.glob(data_path + "/merged_standards/m*.fits")
standard_table = pd.read_csv('./standard_table_v2.txt', index_col=0)  # csv of standards with file and Spectral Type, c/v TBA

# early_k = ('1', '2', '3', '4', '5')
# Only want to look at the K types for the moment
# standard_table_k = standard_table[standard_table['Spectral_Type'].str.startswith(('K1', 'K2', 'K5'))]

standard_shift = standard_table['Wavelength Shift'].values
standard_list = standard_table['Source'].values

hops_table = pd.read_csv('./hops_table.txt')
hops_list = hops_table['Source']

# Determine the maximum length of flux arrays for the standards
max_flux_length = max(len(fits.getdata(file)[1]) for file in standard_list)
max_wavelen_length = max(len(fits.getdata(file)[0]) for file in standard_list)
max_snr_length = max(len(fits.getdata(file)[2]) for file in standard_list)

hops_flux_length = max(len(fits.getdata(file)[1]) for file in hops_table['Source'])
hops_wavelen_length = max(len(fits.getdata(file)[0]) for file in hops_table['Source'])
hops_snr_length = max(len(fits.getdata(file)[2]) for file in hops_table['Source'])

# Initialize stacks with NaN values
wavelen_stack = np.full((max_wavelen_length, len(standard_list)), np.nan)
raw_flux_stack = np.full((max_flux_length, len(standard_list)), np.nan)
snr_stack = np.full((max_snr_length, len(standard_list)), np.nan)
raw_flux_err_stack = np.full((max_snr_length, len(standard_list)), np.nan)

# hops_wavelen_stack = np.full((max_wavelen_length, len(hops_table)), np.nan)
# hops_raw_flux_stack = np.full((max_flux_length, len(hops_table)), np.nan)
# hops_snr_stack = np.full((max_snr_length, len(hops_table)), np.nan)
# hops_raw_flux_err_stack = np.full((max_snr_length, len(hops_table)), np.nan)


# Fill stacks with data
for i, file in enumerate(standard_list):
    # Get wavelength, flux, snr per resolution element data
    data = fits.getdata(file)
    wavelen, flux, snr = data[0], data[1], data[2]
    
    # Clean data a bit
    snr_min = 50  # Minimum SNR
    snr_max = 1e4  # Maximum SNR
    snr_cut = (snr > snr_min)  # & (snr < snr_max)  # Bitwise SNR masking

    flux_min = 10000  # Minimum flux
    flux_cut = flux > flux_min  # Bitwise flux masking

    wavelen_min = 2.0
    wavelen_max = 2.4
    wavelen_cut = (wavelen > wavelen_min) & (wavelen < wavelen_max)

    # Apply masks and remove NaNs and infs
    mask = flux_cut & wavelen_cut & snr_cut
    
    wavelen = wavelen[mask]
    flux = flux[mask]
    snr = snr[mask]

    # Remove NaNs and infs from wavelen, flux, and snr arrays
    valid_indices = ~np.isnan(wavelen) & ~np.isnan(flux) & ~np.isnan(snr)
    
    wavelen = wavelen[valid_indices] - standard_shift[i]
    flux = flux[valid_indices]
    snr = snr[valid_indices]

    # Check for NaNs in the final arrays
    if np.any(np.isnan(wavelen)) or np.any(np.isnan(flux)) or np.any(np.isnan(snr)):
        print(f"NaNs found in data for file {file} after cleaning")

    wavelen_stack[:len(wavelen), i] = wavelen  # Wavelength arrays for each standard
    raw_flux_stack[:len(flux), i] = flux
    snr_stack[:len(snr), i] = snr
    raw_flux_err_stack[:len(flux), i] = flux / snr

# for i, file in enumerate(hops_list):
#     # Get wavelength, flux, snr per resolution element data
#     data = fits.getdata(file)
#     wavelen, flux, snr = data[0], data[1], data[2]
    
#     # Clean data a bit
#     snr_min = 50  # Minimum SNR
#     snr_max = 1e4  # Maximum SNR
#     snr_cut = (snr > snr_min)  # & (snr < snr_max)  # Bitwise SNR masking

#     flux_min = 10000  # Minimum flux
#     flux_cut = flux > flux_min  # Bitwise flux masking

#     wavelen_min = 2.0
#     wavelen_max = 2.4
#     wavelen_cut = (wavelen > wavelen_min) & (wavelen < wavelen_max)

#     # Apply masks and remove NaNs and infs
#     mask = flux_cut & wavelen_cut & snr_cut
    
#     wavelen = wavelen[mask]
#     flux = flux[mask]
#     snr = snr[mask]

#     # Remove NaNs and infs from wavelen, flux, and snr arrays
#     valid_indices = ~np.isnan(wavelen) & ~np.isnan(flux) & ~np.isnan(snr)
    
#     wavelen = wavelen[valid_indices] - standard_shift[i]
#     flux = flux[valid_indices]
#     snr = snr[valid_indices]

#     # Check for NaNs in the final arrays
#     if np.any(np.isnan(wavelen)) or np.any(np.isnan(flux)) or np.any(np.isnan(snr)):
#         print(f"NaNs found in data for file {file} after cleaning")

#     hops_wavelen_stack[:len(wavelen), i] = wavelen  # Wavelength arrays for each standard
#     hops_raw_flux_stack[:len(flux), i] = flux
#     hops_snr_stack[:len(snr), i] = snr
#     hops_raw_flux_err_stack[:len(flux), i] = flux / snr

# Directly query NIST to find line features in K-band
with warnings.catch_warnings():  # Ignore warnings
    warnings.simplefilter('ignore') 
    lines_table = Nist.query(2.08*u.um,2.35*u.um,
                    linename = 'Na I, Sc I, Si I, Fe I, Fe II, Al I, Mg I, Ca I, H I, Ti I',
                    energy_level_unit='eV',output_order='wavelength')

igrins_wav_cut = (lines_table['Observed'] > 2.08) & (lines_table['Observed'] < 2.35)
lines_table = lines_table[igrins_wav_cut]
# lines_table = pd.read_csv('lines_table.txt')

# Make masks for the table of all the lines just in case I want to peek at certain transitions/wavelengths
na1_mask = lines_table['Spectrum'] == 'Na I'
sc1_mask = lines_table['Spectrum'] == 'Sc I'
si1_mask = lines_table['Spectrum'] == 'Si I'
fe1_mask = lines_table['Spectrum'] == 'Fe I'
fe2_mask = lines_table['Spectrum'] == 'Fe II'
al1_mask = lines_table['Spectrum'] == 'Al I'
mg1_mask = lines_table['Spectrum'] == 'Mg I'
ca1_mask = lines_table['Spectrum'] == 'Ca I'
h1_mask  = lines_table['Spectrum'] == 'H I'
ti1_mask = lines_table['Spectrum'] == 'Ti I'

# Just add all the masks to a list for the sake of my plotting a few cells down
mask_list = [na1_mask,sc1_mask,si1_mask,fe1_mask,al1_mask,mg1_mask,ca1_mask,h1_mask,ti1_mask]
color_list = ['purple', 'orange', 'green', 'blue', 'brown', 'crimson', 'olive', 'cyan', 'darkgreen']

Veiling (?) Lines

Fe I ~ 2.084 (lines_table[fe1_mask][6]), 2.275 (lines_table[fe1_mask][109]), 2.284 (lines_table[fe1_mask][115])

Si I ~ 2.092 (lines_table[si1_mask][0])

Mg I ~ 2.106 (lines_table[mg1_mask][0]), 2.282 lines_table[mg1_mask][11]

Al I ~ 2.11 (lines_table[al1_mask][1])

Ca I ~ 2.283

Ti I ~ 2.29 (lines_table[ti1_mask][69])

CO (2-0) > 2.294

lines_table[fe1_mask][6]

Row index=6

Spectrum	Observed	Ritz	Transition	Rel.	Aki	fik	Acc.	Ei Ek	Lower level	Upper level	Type	TP	Line
str5	float64	str18	float64	str17	float64	float64	str3	str40	str45	str59	str5	str17	str20
Fe I	2.08465	2.0846493	4796.968	1200	--	--	--	6.01718438 - 6.61193286	3d7.(4F).5s | e 3F | 3	3d7.(4F).5p | 3F* | 4	--	--	L11631

Fe I 2.08465 $\mu m$

line_name = lines_table[fe1_mask][0]['Spectrum']
line_center = lines_table[fe1_mask][6]['Observed']

# from igrins_mod import local_continuum_fit
continuum_stack = []
norm_flux_stack = []
reg_idx_stack = []

# regions [left point,width]
regions = [(-200,10), (200,10)]

# number of regions I use and subtract 1 since index starts at 0
n = len(regions)-1
for i in range(len(standard_table)):
    fig = plt.figure(figsize=(15,3))
    continuum, region_indices = ig.local_continuum_fit(wavelen_stack[:,i],
                                                    raw_flux_stack[:,i],
                                                    poly_order = 1,
                                                    line_center = line_center,
                                                    spec_res = spec_res,
                                                    regions = regions)
    # append indices to list of indices
    reg_idx_stack.append(region_indices)
    
    # append to list of the local continuum arrays
    continuum_stack.append(continuum)

    # normalize flux by dividing the flux by the continuum and create a list
    norm_flux = raw_flux_stack[region_indices[0][0]:region_indices[n][1],i]/continuum_stack[i]
    norm_flux_stack.append(norm_flux)

    plt.plot(wavelen_stack[region_indices[0][0]:region_indices[n][1],i], raw_flux_stack[region_indices[0][0]:region_indices[n][1],i])
    plt.plot(wavelen_stack[region_indices[0][0]:region_indices[n][1],i], continuum_stack[i])
    for j in range(len(region_indices)):
        plt.axvspan(wavelen_stack[region_indices[n-j][0],i],wavelen_stack[region_indices[n-j][1],i],color='red',alpha=.2)

    plt.title(rf"{line_name} {line_center}$\mu m$ {standard_table['Name'][i]} {standard_table['Spectral_Type'][i]}")

    plt.show()

amplitude1_stack = []
center1_stack = []
sigma_stack = []

pcov_stack = []

best_model_stack = []

flux_constant = np.linspace(0,-1,len(standard_table))

# Define initial parameters for Gaussian fitting
init_params = (0.1, lines_table[fe1_mask][6]['Observed'], spec_res)


# wavelen_stack[region_indices[0][0]:region_indices[2][1],i],raw_flux_stack[region_indices[0][0]:region_indices[2][1],i]

for i in range(len(standard_table)):    
    # Perform Gaussian fitting for the current source
    popt, pcov, param_err, best_model = ig.model_fit(ig.gaussian,
                                        wavelen_stack[reg_idx_stack[i][0][0]:reg_idx_stack[i][n][1],i],
                                        norm_flux_stack[i], raw_flux_err_stack[reg_idx_stack[i][0][0]:reg_idx_stack[i][n][1],i],
                                        init_params,
                                        maxfev = 50000)

    amplitude1_stack.append(popt[0])
    center1_stack.append(popt[1])
    sigma_stack.append(popt[2])

    pcov_stack.append(pcov)
    best_model_stack.append(best_model)

# for i in range(len(standard_table)):
    # fwhm = 2*np.sqrt(2*np.log(2))*sigma1_stack[i]
    # area = np.abs(amplitude_stack[i]*sigma_stack[i]*np.sqrt(2*np.pi))
    # ew = np.trapz(1-best_model_stack[i],wavelen_stack[reg_idx_stack[i][0][0]:reg_idx_stack[i][n][1],i])
    
fe1_2_0846_fit = []

for i in range(len(standard_table)):
    fe1_2_0846_fit.append(ig.gaussian(wavelen_stack[reg_idx_stack[i][0][0]:reg_idx_stack[i][n][1],i],amplitude1_stack[i],center1_stack[i],sigma_stack[i]))
    
    area_int = trapz(1-best_model_stack[i],wavelen_stack[reg_idx_stack[i][0][0]:reg_idx_stack[i][n][1],i])
    area_v2 = np.abs(amplitude1_stack[i])*sigma_stack[i]*np.sqrt(2*np.pi)

    # print(f'{area_int} {area_v2}')

ew1_stack = [] # empty list to load in equivalent widths

# Integate 1-best_model to get the areas (equivalent widths) of each component Gaussian
for i in range(len(standard_table)):
    fig = plt.figure(figsize=(15,3))

    plt.plot(wavelen_stack[reg_idx_stack[i][0][0]:reg_idx_stack[i][n][1],i],norm_flux_stack[i], c='black')
    plt.plot(wavelen_stack[reg_idx_stack[i][0][0]:reg_idx_stack[i][n][1],i],best_model_stack[i],c='red',lw=2,label='Best Model')
    

    # plt.fill_between(wavelen_stack[reg_idx_stack[i][0][0]:reg_idx_stack[i][n][1],i],best_model_stack[i],0,color='red',alpha=0.2)
    # Fe I
    # plt.plot(wavelen_stack[reg_idx_stack[i][0][0]:reg_idx_stack[i][n][1],i],fe1_2_0846_fit[i])#, label = 'Fe I 2.2263 Fit')
    # plt.ylim(0)
    plt.ylabel('Flux + Constant')
    plt.xlabel(r'Wavelength [$\mu$m]')

    plt.title(f"{line_name} ({line_center}$\mu$m) {standard_table['Name'][i]} {standard_table['Spectral_Type'][i]}")

    plt.legend(loc='best')
    plt.show()

# Define the region for fitting
line_name = lines_table[fe1_mask][0]['Spectrum'] # Species
line_center = lines_table[fe1_mask][109]['Observed'] # Wavelength

# from igrins_mod import local_continuum_fit
continuum_stack = []
norm_flux_stack = []
reg_idx_stack = []

# regions [left point,width]
regions = [(-200,10), (200,10)]

# number of regions I use and subtract 1 since index starts at 0
n = len(regions)-1
for i in range(len(standard_table)):
    fig = plt.figure(figsize=(15,3))
    continuum, region_indices = ig.local_continuum_fit(wavelen_stack[:,i],
                                                    raw_flux_stack[:,i],
                                                    poly_order = 1,
                                                    line_center = line_center,
                                                    spec_res = spec_res,
                                                    regions = regions)
    # append indices to list of indices
    reg_idx_stack.append(region_indices)
    
    # append to list of the local continuum arrays
    continuum_stack.append(continuum)

    # normalize flux by dividing the flux by the continuum and create a list
    norm_flux = raw_flux_stack[region_indices[0][0]:region_indices[n][1],i]/continuum_stack[i]
    norm_flux_stack.append(norm_flux)


    plt.plot(wavelen_stack[region_indices[0][0]:region_indices[n][1],i], raw_flux_stack[region_indices[0][0]:region_indices[n][1],i])
    plt.plot(wavelen_stack[region_indices[0][0]:region_indices[n][1],i], continuum_stack[i])
    for j in range(len(region_indices)):
        plt.axvspan(wavelen_stack[region_indices[n-j][0],i],wavelen_stack[region_indices[n-j][1],i],color='red',alpha=.2)

    plt.title(f"{standard_table['Name'][i]} {standard_table['Spectral_Type'][i]}")

    plt.show()

amplitude1_stack = []
center1_stack = []
sigma_stack = []

pcov_stack = []

best_model_stack = []

flux_constant = np.linspace(0,-1,len(standard_list))

# Define initial parameters for Gaussian fitting
init_params = (0.1, lines_table[fe1_mask][6]['Observed'], spec_res)


# wavelen_stack[region_indices[0][0]:region_indices[2][1],i],raw_flux_stack[region_indices[0][0]:region_indices[2][1],i]

for i in range(len(standard_table)):    
    # Perform Gaussian fitting for the current source
    popt, pcov, param_err, best_model = ig.model_fit(ig.gaussian,wavelen_stack[reg_idx_stack[i][0][0]:reg_idx_stack[i][n][1],i],
                                        norm_flux_stack[i], raw_flux_err_stack[reg_idx_stack[i][0][0]:reg_idx_stack[i][n][1],i],
                                        init_params,
                                        maxfev=5000)

    amplitude1_stack.append(popt[0])
    center1_stack.append(popt[1])
    sigma_stack.append(popt[2])

    pcov_stack.append(pcov)
    best_model_stack.append(best_model)

# for i in range(len(standard_table)):
    # fwhm = 2*np.sqrt(2*np.log(2))*sigma1_stack[i]
    # area = np.abs(amplitude_stack[i]*sigma_stack[i]*np.sqrt(2*np.pi))
    # ew = np.trapz(1-best_model_stack[i],wavelen_stack[reg_idx_stack[i][0][0]:reg_idx_stack[i][n][1],i])
    
# fe1_2_0846_fit = []

# for i in range(len(standard_table)):
#     fe1_2_0846_fit.append(ig.gaussian(wavelen_stack[reg_idx_stack[i][0][0]:reg_idx_stack[i][n][1],i],amplitude1_stack[i],center1_stack[i],sigma_stack[i]))
    
#     area_int = trapz(1-best_model_stack[i],wavelen_stack[reg_idx_stack[i][0][0]:reg_idx_stack[i][n][1],i])
#     area_v2 = np.abs(amplitude1_stack[i])*sigma_stack[i]*np.sqrt(2*np.pi)

#     # print(f'{area_int} {area_v2}')

# ew1_stack = [] # empty list to load in equivalent widths

# Integate 1-best_model to get the areas (equivalent widths) of each component Gaussian
for i in range(len(standard_table)):
    fig = plt.figure(figsize=(15,3))

    plt.plot(wavelen_stack[reg_idx_stack[i][0][0]:reg_idx_stack[i][n][1],i],norm_flux_stack[i], c='black')
    plt.plot(wavelen_stack[reg_idx_stack[i][0][0]:reg_idx_stack[i][n][1],i],best_model_stack[i],c='red',lw=2,label='Best Model')
    

    # plt.fill_between(wavelen_stack[reg_idx_stack[i][0][0]:reg_idx_stack[i][n][1],i],best_model_stack[i],0,color='red',alpha=0.2)
    # Fe I
    # plt.plot(wavelen_stack[reg_idx_stack[i][0][0]:reg_idx_stack[i][n][1],i],fe1_2_0846_fit[i])#, label = 'Fe I 2.2263 Fit')
    # plt.ylim(0)
    plt.ylabel('Flux + Constant')
    plt.xlabel(r'Wavelength [$\mu$m]')

    plt.title(f"{line_name} ({line_center}$\mu$m) {standard_table['Name'][i]} {standard_table['Spectral_Type'][i]}")

    plt.legend(loc='best')
    plt.show()

c:\Users\Savio\anaconda3\envs\muler_dev\lib\site-packages\scipy\optimize\_minpack_py.py:1010: OptimizeWarning: Covariance of the parameters could not be estimated
  warnings.warn('Covariance of the parameters could not be estimated',

Si 1

# Define the region for fitting
line_name = lines_table[si1_mask][0]['Spectrum'] # Species
line_center = lines_table[si1_mask][0]['Observed'] # Wavelength

# from igrins_mod import local_continuum_fit
continuum_stack = []
norm_flux_stack = []
reg_idx_stack = []

# regions [left point,width]
regions = [(-200,10), (200,10)]

# number of regions I use and subtract 1 since index starts at 0
n = len(regions)-1
for i in range(len(standard_table)):
    fig = plt.figure(figsize=(15,3))
    continuum, region_indices = ig.local_continuum_fit(wavelen_stack[:,i],
                                                    raw_flux_stack[:,i],
                                                    poly_order = 1,
                                                    line_center = line_center,
                                                    spec_res = spec_res,
                                                    regions = regions)
    # append indices to list of indices
    reg_idx_stack.append(region_indices)
    
    # append to list of the local continuum arrays
    continuum_stack.append(continuum)

    # normalize flux by dividing the flux by the continuum and create a list
    norm_flux = raw_flux_stack[region_indices[0][0]:region_indices[n][1],i]/continuum_stack[i]
    norm_flux_stack.append(norm_flux)


    plt.plot(wavelen_stack[region_indices[0][0]:region_indices[n][1],i], raw_flux_stack[region_indices[0][0]:region_indices[n][1],i])
    plt.plot(wavelen_stack[region_indices[0][0]:region_indices[n][1],i], continuum_stack[i])
    for j in range(len(region_indices)):
        plt.axvspan(wavelen_stack[region_indices[n-j][0],i],wavelen_stack[region_indices[n-j][1],i],color='red',alpha=.2)

    plt.title(f"{line_name} {line_center} {standard_table['Name'][i]} {standard_table['Spectral_Type'][i]}")

    plt.show()

Mg I

# Define the region for fitting
line_name = lines_table[mg1_mask][0]['Spectrum'] # Species
line_center = lines_table[mg1_mask][0]['Observed'] # Wavelength

# from igrins_mod import local_continuum_fit
continuum_stack = []
norm_flux_stack = []
reg_idx_stack = []

# regions [left point,width]
regions = [(-200,10), (200,10)]

# number of regions I use and subtract 1 since index starts at 0
n = len(regions)-1
for i in range(len(standard_table)):
    fig = plt.figure(figsize=(15,3))
    continuum, region_indices = ig.local_continuum_fit(wavelen_stack[:,i],
                                                    raw_flux_stack[:,i],
                                                    poly_order = 1,
                                                    line_center = line_center,
                                                    spec_res = spec_res,
                                                    regions = regions)
    # append indices to list of indices
    reg_idx_stack.append(region_indices)
    
    # append to list of the local continuum arrays
    continuum_stack.append(continuum)

    # normalize flux by dividing the flux by the continuum and create a list
    norm_flux = raw_flux_stack[region_indices[0][0]:region_indices[n][1],i]/continuum_stack[i]
    norm_flux_stack.append(norm_flux)


    plt.plot(wavelen_stack[region_indices[0][0]:region_indices[n][1],i], raw_flux_stack[region_indices[0][0]:region_indices[n][1],i])
    plt.plot(wavelen_stack[region_indices[0][0]:region_indices[n][1],i], continuum_stack[i])
    for j in range(len(region_indices)):
        plt.axvspan(wavelen_stack[region_indices[n-j][0],i],wavelen_stack[region_indices[n-j][1],i],color='red',alpha=.2)

    plt.title(f"{standard_table['Name'][i]} {standard_table['Spectral_Type'][i]}")

    plt.show()

Al I

# Define the region for fitting
line_name = lines_table[al1_mask][0]['Spectrum'] # Species
line_center = lines_table[al1_mask][1]['Observed'] # Wavelength

# from igrins_mod import local_continuum_fit
continuum_stack = []
norm_flux_stack = []
reg_idx_stack = []

# regions [left point,width]
regions = [(-200,10), (200,10)]

# number of regions I use and subtract 1 since index starts at 0
n = len(regions)-1
for i in range(len(standard_table)):
    fig = plt.figure(figsize=(15,3))
    continuum, region_indices = ig.local_continuum_fit(wavelen_stack[:,i],
                                                    raw_flux_stack[:,i],
                                                    poly_order = 1,
                                                    line_center = line_center,
                                                    spec_res = spec_res,
                                                    regions = regions)
    # append indices to list of indices
    reg_idx_stack.append(region_indices)
    
    # append to list of the local continuum arrays
    continuum_stack.append(continuum)

    # normalize flux by dividing the flux by the continuum and create a list
    norm_flux = raw_flux_stack[region_indices[0][0]:region_indices[n][1],i]/continuum_stack[i]
    norm_flux_stack.append(norm_flux)


    plt.plot(wavelen_stack[region_indices[0][0]:region_indices[n][1],i], raw_flux_stack[region_indices[0][0]:region_indices[n][1],i])
    plt.plot(wavelen_stack[region_indices[0][0]:region_indices[n][1],i], continuum_stack[i])
    for j in range(len(region_indices)):
        plt.axvspan(wavelen_stack[region_indices[n-j][0],i],wavelen_stack[region_indices[n-j][1],i],color='red',alpha=.2)

    plt.title(f"{standard_table['Name'][i]} {standard_table['Spectral_Type'][i]}")

    plt.show()

Ti I

# Define the region for fitting
line_name = lines_table[ti1_mask][0]['Spectrum'] # Species
line_center = lines_table[ti1_mask][69]['Observed'] # Wavelength

# from igrins_mod import local_continuum_fit
continuum_stack = []
norm_flux_stack = []
reg_idx_stack = []

# regions [left point,width]
regions = [(-200,10), (200,10)]

# number of regions I use and subtract 1 since index starts at 0
n = len(regions)-1
for i in range(len(standard_table)):
    fig = plt.figure(figsize=(15,3))
    continuum, region_indices = ig.local_continuum_fit(wavelen_stack[:,i],
                                                    raw_flux_stack[:,i],
                                                    poly_order = 1,
                                                    line_center = line_center,
                                                    spec_res = spec_res,
                                                    regions = regions)
    # append indices to list of indices
    reg_idx_stack.append(region_indices)
    
    # append to list of the local continuum arrays
    continuum_stack.append(continuum)

    # normalize flux by dividing the flux by the continuum and create a list
    norm_flux = raw_flux_stack[region_indices[0][0]:region_indices[n][1],i]/continuum_stack[i]
    norm_flux_stack.append(norm_flux)


    plt.plot(wavelen_stack[region_indices[0][0]:region_indices[n][1],i], raw_flux_stack[region_indices[0][0]:region_indices[n][1],i])
    plt.plot(wavelen_stack[region_indices[0][0]:region_indices[n][1],i], continuum_stack[i])
    for j in range(len(region_indices)):
        plt.axvspan(wavelen_stack[region_indices[n-j][0],i],wavelen_stack[region_indices[n-j][1],i],color='red',alpha=.2)

    plt.title(f"{standard_table['Name'][i]} {standard_table['Spectral_Type'][i]}")

    plt.show()

Temperature Sensitive Lines

Sc I (2.2058 $\mu m$ and 2.2071 $\mu m$): index 1 and 2

Si I (2.2068 $\mu m$)

Mg/Al (~2.11 $\mu m$)

Fe I (2.2205–2.2346 $\mu m$) Over Ti region: index 10 and 11

Ca I (2.2614, 2.2631, 2.2657 $\mu m$): Only three lines

Na Interval ~2.21

Na I Doublet (2.2062, 2.2089) & Sc I (2.2058)

# Define the region for fitting
line_name = lines_table[na1_mask][0]['Spectrum'] # Species
line_center = lines_table[na1_mask][0]['Observed'] # Wavelength

# from igrins_mod import local_continuum_fit
continuum_stack = []
norm_flux_stack = []
reg_idx_stack = []

# regions [left point,width]
regions = [(-200,30),(-140,30),(150,40),(350,30)]

# number of regions I use and subtract 1 since index starts at 0
n = len(regions)-1
poly_deg = 3
for i in range(len(standard_table)):
    fig = plt.figure(figsize=(15,3))
    continuum, region_indices = ig.local_continuum_fit(wavelen_stack[:,i],
                                                    raw_flux_stack[:,i],
                                                    poly_order = poly_deg,
                                                    line_center = line_center,
                                                    spec_res = spec_res,
                                                    regions = regions)
    # append indices to list of indices
    reg_idx_stack.append(region_indices)
    
    # append to list of the local continuum arrays
    continuum_stack.append(continuum)

    # normalize flux by dividing the flux by the continuum and create a list
    norm_flux = raw_flux_stack[region_indices[0][0]:region_indices[n][1],i]/continuum_stack[i]
    norm_flux_stack.append(norm_flux)


    plt.plot(wavelen_stack[region_indices[0][0]:region_indices[n][1],i], raw_flux_stack[region_indices[0][0]:region_indices[n][1],i])
    plt.plot(wavelen_stack[region_indices[0][0]:region_indices[n][1],i], continuum_stack[i])

    for j in range(len(region_indices)):
        plt.axvspan(wavelen_stack[region_indices[n-j][0],i],wavelen_stack[region_indices[n-j][1],i],color='red',alpha=.2)

    plt.title(f"Na Interval {standard_table['Name'][i]} {standard_table['Spectral_Type'][i]}")

    plt.show()

# Initialize storage lists
params_stack = [] # list of best of fit parameters
params_error_stack = [] # list of errors on each parameter
best_model_stack = [] # list of best fits to the dats
result_stack = [] # lmfit results
dely_stack = [] # lmfit errors of the fits

for i in range(len(standard_table)):
    # Define initial parameters for Gaussian fitting
    params = Parameters()
    params.add('amp1', value=-1e-5, min=-1, max=0)
    params.add('c1', value=lines_table[sc1_mask][19]['Observed'], min=lines_table[sc1_mask][19]['Observed']-0.5e-4,max=lines_table[sc1_mask][19]['Observed']+0.5e-4)
    params.add('std1', value=spec_res)

    params.add('amp2', value=-1e-5, min=-1, max=0)
    params.add('c2', value=line_center, min=line_center-0.5e-4,max=line_center+0.5e-4)
    params.add('std2', value=spec_res)

    params.add('amp3', value=-1e-5, min=-1, max=0)
    params.add('c3', value=lines_table[si1_mask][2]['Observed'],min=lines_table[si1_mask][2]['Observed']-0.5e-4,max=lines_table[si1_mask][2]['Observed']+0.5e-4)
    params.add('std3', value=spec_res)

    params.add('amp4', value=-1e-5, min=-2e-5, max=0)
    params.add('c4', value=lines_table[sc1_mask][20]['Observed'],min=lines_table[sc1_mask][20]['Observed']-0.5e-4,max=lines_table[sc1_mask][20]['Observed']+0.5e-4)
    params.add('std4', value=spec_res)

    params.add('amp5', value=-1e-5, min=-1, max=0)
    params.add('c5', value=lines_table[na1_mask][1]['Observed'],min=lines_table[na1_mask][1]['Observed']-0.5e-4,max=lines_table[na1_mask][1]['Observed']+0.5e-4)
    params.add('std5', value=spec_res)


    # mathematical/logic constraints on parameters
    # enforce same widths for the Scandium and Sodium lines
    # params['std4'].expr = 'std1'
    # params['std5'].expr = 'std2'

    # Adding a parameter for the wavelength differences
    # params.add('delta', value=1e-3)

    # Constraining the center parameters to have the same difference
    # params['c2'].expr = 'c1 + delta'
    # params['c3'].expr = 'c1 + 2*delta'
    # params['c4'].expr = 'c1 + 3*delta'
    # params['c5'].expr = 'c1 + 4*delta'

    # Define the model
    model = Model(ig.five_gaussian, nan_policy='omit')

    # Perform the fit
    result = model.fit(norm_flux_stack[i], params,
                       x=wavelen_stack[reg_idx_stack[i][0][0]:reg_idx_stack[i][n][1], i],
                       weights=1 / raw_flux_err_stack[reg_idx_stack[i][0][0]:reg_idx_stack[i][n][1], i])

    result_stack.append(result) # lmfit ModelResult objects
    dely = result.eval_uncertainty(sigma=5) # lmfit 
    dely_stack.append(dely)

    # Append the best_model, error and parameters stacks
    # could probably make a loop to the tune of for value in result etc...
    params_stack.extend([result.params['amp1'].value, result.params['c1'].value, result.params['std1'].value,
                         result.params['amp2'].value, result.params['c2'].value, result.params['std2'].value,
                         result.params['amp3'].value, result.params['c3'].value, result.params['std3'].value,
                         result.params['amp4'].value, result.params['c4'].value, result.params['std4'].value,
                         result.params['amp5'].value, result.params['c5'].value, result.params['std5'].value])
    
    params_error_stack.extend([result.params['amp1'].stderr, result.params['c1'].stderr, result.params['std1'].stderr,
                               result.params['amp2'].stderr, result.params['c2'].stderr, result.params['std2'].stderr,
                               result.params['amp3'].stderr, result.params['c3'].stderr, result.params['std3'].stderr,
                               result.params['amp4'].stderr, result.params['c4'].stderr, result.params['std4'].stderr,
                               result.params['amp5'].stderr, result.params['c5'].stderr, result.params['std5'].stderr])

    best_model_stack.append(result.best_fit)
    # print(result.fit_report())

params_arr = np.array(params_stack)

amps = params_arr[0::3]
amps_err = np.array(params_error_stack)[0::3]

centers = params_arr[1::3]
center_err = np.array(params_error_stack)[1::3]

sigmas = params_arr[2::3]
sigmas_err = np.array(params_error_stack)[2::3]

# Generate the fits for each Gaussian component
na1_2_2062_fit = []
sc1_2_2058_fit = []
si1_2_2068_fit = []
sc1_2_2071_fit = []
na1_2_2089_fit = []
num_gauss = 5

for i in range(len(standard_table)):
    sc1_2_2058_fit.append(ig.gaussian(wavelen_stack[reg_idx_stack[i][0][0]:reg_idx_stack[i][n][1], i], amps[0::num_gauss][i], centers[0::num_gauss][i], sigmas[0::num_gauss][i]))
    na1_2_2062_fit.append(ig.gaussian(wavelen_stack[reg_idx_stack[i][0][0]:reg_idx_stack[i][n][1], i], amps[1::num_gauss][i], centers[1::num_gauss][i], sigmas[1::num_gauss][i]))
    si1_2_2068_fit.append(ig.gaussian(wavelen_stack[reg_idx_stack[i][0][0]:reg_idx_stack[i][n][1], i], amps[2::num_gauss][i], centers[2::num_gauss][i], sigmas[2::num_gauss][i]))
    sc1_2_2071_fit.append(ig.gaussian(wavelen_stack[reg_idx_stack[i][0][0]:reg_idx_stack[i][n][1], i], amps[3::num_gauss][i], centers[3::num_gauss][i], sigmas[3::num_gauss][i]))
    na1_2_2089_fit.append(ig.gaussian(wavelen_stack[reg_idx_stack[i][0][0]:reg_idx_stack[i][n][1], i], amps[4::num_gauss][i], centers[4::num_gauss][i], sigmas[4::num_gauss][i]))

for i in range(len(standard_table)):
    # Create subplots with adjusted size ratios
    fig, (ax1, ax2) = plt.subplots(2, 1, figsize=(15, 6), gridspec_kw={'height_ratios': [4, 1]}, sharex=True)

    ax1.plot(wavelen_stack[reg_idx_stack[i][0][0]:reg_idx_stack[i][n][1], i], norm_flux_stack[i], c='black', alpha=0.2, label='Data')
    ax1.plot(wavelen_stack[reg_idx_stack[i][0][0]:reg_idx_stack[i][n][1], i], best_model_stack[i], c='red', lw=1, label='Best Model')
    # ax1.errorbar(wavelen_stack[reg_idx_stack[i][0][0]:reg_idx_stack[i][n][1], i], best_model_stack[i], yerr=dely_stack[i],fmt='ro', ms=3, lw=1, label='Best Model')
    
    ax1.fill_between(wavelen_stack[reg_idx_stack[i][0][0]:reg_idx_stack[i][n][1], i],
        best_model_stack[i]-dely_stack[i], best_model_stack[i]+dely_stack[i],
        alpha=0.2, color='red')

    # plot the locations of fitted centers
    ax1.scatter(y=[1.03] * len(centers[i * num_gauss:(i + 1) * num_gauss]), x=centers[i * num_gauss:(i + 1) * num_gauss], marker="|", s=100, color='black')

    ax1.plot(wavelen_stack[reg_idx_stack[i][0][0]:reg_idx_stack[i][n][1], i], sc1_2_2058_fit[i], ls='-.',
             label=r'$\bf{Sc~I~2.2058}$' "\n" rf"A={amps[0::num_gauss][i]:.3e}" "\n" rf"$\mu$={centers[0::num_gauss][i]:.5f}" "\n" rf"$\sigma$= {sigmas[0::num_gauss][i]:.3e}")

    ax1.plot(wavelen_stack[reg_idx_stack[i][0][0]:reg_idx_stack[i][n][1], i], na1_2_2062_fit[i], ls='-.',
             label=r'$\bf{Na~I~2.2062}$' "\n" rf"A={amps[1::num_gauss][i]:.3e}" "\n" rf"$\mu$={centers[1::num_gauss][i]:.5f}" "\n" rf"$\sigma$= {sigmas[1::num_gauss][i]:.3e}")

    ax1.plot(wavelen_stack[reg_idx_stack[i][0][0]:reg_idx_stack[i][n][1], i], si1_2_2068_fit[i], ls='-.',
             label=r'$\bf{Si~I~2.2068}$' "\n" rf"A={amps[2::num_gauss][i]:.3e}" "\n" rf"$\mu$={centers[2::num_gauss][i]:.5f}" "\n" rf"$\sigma$= {sigmas[2::num_gauss][i]:.3e}")
            
    ax1.plot(wavelen_stack[reg_idx_stack[i][0][0]:reg_idx_stack[i][n][1], i], sc1_2_2071_fit[i], ls='-.',
             label=r'$\bf{Sc~I~2.2071}$' "\n" rf"A={amps[3::num_gauss][i]:.3e}" "\n" rf"$\mu$={centers[3::num_gauss][i]:.5f}" "\n" rf"$\sigma$= {sigmas[3::num_gauss][i]:.3e}")

    ax1.plot(wavelen_stack[reg_idx_stack[i][0][0]:reg_idx_stack[i][n][1], i], na1_2_2089_fit[i], ls='-.',
             label=r'$\bf{Na~I~2.2089}$' "\n" rf"A={amps[4::num_gauss][i]:.3e}" "\n" rf"$\mu$={centers[4::num_gauss][i]:.5f}" "\n" rf"$\sigma$= {sigmas[4::num_gauss][i]:.3e}")
    # ax1.set_ylim(top=1.05)
    ax1.set_ylabel('Flux + Constant')

    # Calculate residuals
    residuals = norm_flux_stack[i] - best_model_stack[i]
    sqsum_res = np.sum(residuals**2)
    sqsum = np.sum((norm_flux_stack[i] - np.mean(norm_flux_stack[i]))**2)
    R2 = 1 - (sqsum_res / sqsum)

    # Plot residuals
    ax2.scatter(wavelen_stack[reg_idx_stack[i][0][0]:reg_idx_stack[i][n][1], i], residuals, c='red', label='Residuals', s=3)
    ax2.axhline(y=0, color='grey', linestyle='--')
    ax2.set_xlabel(r'Wavelength [$\mu$m]')
    ax2.set_ylabel('res')
    ax2.set_title(rf'$R^2$={R2:.4f}', fontsize=12)

    plt.suptitle(rf"Na Interval {standard_table['Name'][i]} {standard_table['Spectral_Type'][i]}")
    ax1.legend(bbox_to_anchor=(1, 1), fontsize=12)
    plt.show()
    # print(result_stack[i].fit_report())

# params_stack = []
# params_error_stack = []

# pcov_stack = []
# best_model_stack = []

# for i in range(len(standard_table)):
#     # Define initial parameters for Gaussian fitting
#     init_params = (-0.05, lines_table[sc1_mask][19]['Observed'], spec_res,
#                           -0.5, line_center, spec_res,
#                           -0.05, lines_table[si1_mask][2]['Observed'], spec_res,
#                           -0.5, lines_table[na1_mask][1]['Observed'], spec_res)

#     # bounds for each parameter
#     # gauss_bounds = ([-1,-np.inf,-np.inf, -1,-np.inf,-np.inf, -1,-np.inf,-np.inf, -1,-np.inf,-np.inf],
#     #                 [0,np.inf,np.inf, 0,np.inf,np.inf, 0,np.inf,np.inf, 0,np.inf,np.inf])

#     # Perform Gaussian fitting for the current source
#     # add error keyword then propogate for equiv width
#     # compare equiv widths between objects of the same spectral types
#     popt, pcov, param_err, best_model = ig.model_fit(ig.four_gaussian,
#                                         wavelen_stack[reg_idx_stack[i][0][0]:reg_idx_stack[i][n][1],i],
#                                         norm_flux_stack[i],
#                                         raw_flux_err_stack[reg_idx_stack[i][0][0]:reg_idx_stack[i][n][1],i],
#                                         init_params,
#                                         maxfev=5000)

#     # Append the best_model, error and parameters stacks
#     params_error_stack.extend(param_err)

#     params_stack.extend(popt)
#     params_arr = np.array(params_stack)

#     # best fit amplitudes
#     amps = params_arr[0::3]
#     amps_err = params_error_stack[0::3]

#     # best fit centers(means)
#     centers = params_arr[1::3]
#     center_err = params_error_stack[1::3]

#     # best fit widths
#     sigmas = params_arr[2::3]
#     sigmas_err = params_error_stack[2::3]

#     pcov_stack.append(pcov)
#     best_model_stack.append(best_model)

# na1_2_2062_fit = []
# sc1_2_2058_fit = []
# si1_2_2068_fit = []
# na1_2_2089_fit = []

# for i in range(len(standard_table)):
#     sc1_2_2058_fit.append(ig.gaussian(wavelen_stack[reg_idx_stack[i][0][0]:reg_idx_stack[i][n][1],i],amps[0::4][i],centers[0::4][i],sigmas[0::4][i]))
#     na1_2_2062_fit.append(ig.gaussian(wavelen_stack[reg_idx_stack[i][0][0]:reg_idx_stack[i][n][1],i],amps[1::4][i],centers[1::4][i],sigmas[1::4][i]))
#     si1_2_2068_fit.append(ig.gaussian(wavelen_stack[reg_idx_stack[i][0][0]:reg_idx_stack[i][n][1],i],amps[2::4][i],centers[2::4][i],sigmas[2::4][i]))
#     na1_2_2089_fit.append(ig.gaussian(wavelen_stack[reg_idx_stack[i][0][0]:reg_idx_stack[i][n][1],i],amps[3::4][i],centers[3::4][i],sigmas[3::4][i]))

# for i in range(len(standard_table)):
#     # Create subplots with adjusted size ratios
#     fig, (ax1, ax2) = plt.subplots(2, 1, figsize=(15, 6), gridspec_kw={'height_ratios': [4, 1]}, sharex=True)

#     ax1.plot(wavelen_stack[reg_idx_stack[i][0][0]:reg_idx_stack[i][n][1],i], norm_flux_stack[i], c='black',alpha=0.2, label='Data')
    
#     ax1.scatter(wavelen_stack[reg_idx_stack[i][0][0]:reg_idx_stack[i][n][1], i], best_model_stack[i], color='red', label='Best Model')
#     ax1.scatter(y=[1.03] * len(centers[i * 4:(i + 1) * 4]), x=centers[i * 4:(i + 1) * 4], marker="|", s=50, color='black')

#     # for j in range(len(amps[::3])):
#     ax1.plot(wavelen_stack[reg_idx_stack[i][0][0]:reg_idx_stack[i][n][1],i], sc1_2_2058_fit[i],ls='-.',
#         label= r'Sc I 2.2058 Fit' "\n" rf"A={amps[0::4][i]:.3e}" "\n" rf"$\mu$={centers[0::4][i]:.5f}" "\n" rf"$\sigma$= {sigmas[0::4][i]:.3e}")

#     ax1.plot(wavelen_stack[reg_idx_stack[i][0][0]:reg_idx_stack[i][n][1],i], na1_2_2062_fit[i],ls='-.',
#         label= r'Na I 2.2062 Fit' "\n" rf"A={amps[1::4][i]:.3e}" "\n" rf"$\mu$={centers[1::4][i]:.5f}" "\n" rf"$\sigma$= {sigmas[1::4][i]:.3e}")

#     ax1.plot(wavelen_stack[reg_idx_stack[i][0][0]:reg_idx_stack[i][n][1],i], si1_2_2068_fit[i],ls='-.',
#         label = r'Si I 2.2068 Fit' "\n" rf"A={amps[2::4][i]:.3e}" "\n" rf"$\mu$={centers[2::4][i]:.5f}" "\n" rf"$\sigma$= {sigmas[2::4][i]:.3e}")

#     ax1.plot(wavelen_stack[reg_idx_stack[i][0][0]:reg_idx_stack[i][n][1],i], na1_2_2089_fit[i],ls='-.',
#         label= r'Na I 2.2089 Fit' "\n" rf"A={amps[3::4][i]:.3e}" "\n" rf"$\mu$={centers[3::4][i]:.5f}" "\n" rf"$\sigma$= {sigmas[3::4][i]:.3e}")

#     ax1.set_ylim(top=1.05)
#     ax1.set_ylabel('Flux + Constant')
    
#     # Calculate residuals and R^2
#     residuals = (norm_flux_stack[i] - best_model_stack[i])
#     sqsum_res = np.sum(residuals**2)
#     sqsum = np.sum((norm_flux_stack[i]-np.mean(norm_flux_stack[i])**2))
#     R2 = 1-(sqsum_res/sqsum)

#     # Plot residuals
#     ax2.scatter(wavelen_stack[reg_idx_stack[i][0][0]:reg_idx_stack[i][n][1], i], residuals, c='red', label='Residuals', s=3)
#     ax2.axhline(y=0, color='grey', linestyle='--')
#     # ax2.set_ylim(-0.1,0.1)
#     ax2.set_xlabel(r'Wavelength [$\mu$m]')
#     ax2.set_ylabel('res')
#     ax2.set_title(rf'$R^2$={R2:.3f}', fontsize=12)

#     plt.suptitle(rf"{line_name} ({line_center}$\mu$m) {standard_table['Name'][i]} {standard_table['Spectral_Type'][i]}")

#     ax1.legend(bbox_to_anchor=(1,1), fontsize=12)
#     plt.show()

ew1_stack = [] # empty list to load in equivalent widths
ew2_stack = []
ew3_stack = []
ew4_stack = []
ew5_stack = []

# Integate 1-best_model to get the areas (equivalent widths) of each component Gaussian
for i in range(len(standard_table)):

    # ew1 = (ig.gaussian_area(amps[0::5][i],sigmas[0::5][i])) # integrate the gaussian
    # ew1_stack.append(ew1)

    # ew2 = (ig.gaussian_area(amps[1::5][i],sigmas[1::5][i]))
    # ew2_stack.append(ew2)

    # ew3 = (ig.gaussian_area(amps[2::5][i],sigmas[2::5][i]))
    # ew3_stack.append(ew3)

    # ew4 = (ig.gaussian_area(amps[3::5][i],sigmas[3::5][i]))
    # ew4_stack.append(ew4)

    # ew5 = (ig.gaussian_area(amps[4::5][i],sigmas[4::5][i]))
    # ew5_stack.append(ew5)
    ew1 = np.trapz(1-na1_2_2062_fit[i],wavelen_stack[reg_idx_stack[i][0][0]:reg_idx_stack[i][n][1],i]) # integrate the gaussian
    ew1_stack.append(ew1)

    ew2 = np.trapz(1-sc1_2_2058_fit[i],wavelen_stack[reg_idx_stack[i][0][0]:reg_idx_stack[i][n][1],i])
    ew2_stack.append(ew2)

    ew3 = np.trapz(1-si1_2_2068_fit[i],wavelen_stack[reg_idx_stack[i][0][0]:reg_idx_stack[i][n][1],i])
    ew3_stack.append(ew3)

    ew4 = np.trapz(1-sc1_2_2071_fit[i],wavelen_stack[reg_idx_stack[i][0][0]:reg_idx_stack[i][n][1],i])
    ew4_stack.append(ew4)

    ew5 = np.trapz(1-na1_2_2089_fit[i],wavelen_stack[reg_idx_stack[i][0][0]:reg_idx_stack[i][n][1],i])
    ew5_stack.append(ew5)

standard_table.loc[:,'ew_na1_2_2062'] = ew1_stack
standard_table.loc[:,'ew_sc1_2_2058'] = ew2_stack
standard_table.loc[:,'ew_si1_2_2068'] = ew3_stack
standard_table.loc[:,'ew_sc1_2_2071'] = ew4_stack
standard_table.loc[:,'ew_na1_2_2089'] = ew5_stack

fig = plt.figure(figsize=(15,5))

plt.plot(standard_table['Name'], standard_table['ew_sc1_2_2058'],ls='-',marker='o', label=f"Sc I 2.2058 $\mu$m")
plt.plot(standard_table['Name'], standard_table['ew_na1_2_2062'],ls='-',marker='o', label=f"Na I 2.2062 $\mu$m")
plt.plot(standard_table['Name'], standard_table['ew_si1_2_2068'],ls='-',marker='o', label=f"Si I 2.2068 $\mu$m")
plt.plot(standard_table['Name'], standard_table['ew_sc1_2_2071'],ls='-',marker='o', label=f"Sc I 2.2071 $\mu$m")
plt.plot(standard_table['Name'], standard_table['ew_na1_2_2089'],ls='-',marker='o', label=f"Na I 2.2089 $\mu$m")

# plt.scatter(standard_table['Spectral_Type'], area_stack)
# plt.ylim(np.nanmedian(ew_stack)+-0.001,np.nanmedian(ew_stack)+0.001)

# plt.ylim(1e-6,1e-3)
plt.xticks(rotation=45)
plt.xlabel('Spectral Type')
plt.ylabel(r'Equivalent Width [$\mu$m]')
# plt.title(f"{line_name} {line_center} $\mu$m")
# plt.yscale('log')

plt.legend(bbox_to_anchor=(1,1))
plt.show()

Mg/Al Interval ~ 2.11

Mg 2.1065 = lines_table[mg1_mask][0]

Al 2.1098 = lines_table[al1_mask][0]

# Define the region for fitting
line_name = lines_table[mg1_mask][0]['Spectrum'] # Species
line_center = lines_table[mg1_mask][0]['Observed'] # Wavelength

# from igrins_mod import local_continuum_fit
continuum_stack = []
norm_flux_stack = []
reg_idx_stack = []

# regions [left point,width]
regions = [(-160,30),(-80,20),(100,20),(200,20),(270,20),(430,20)]

# number of regions I use and subtract 1 since index starts at 0
n = len(regions)-1
poly_deg = 3
for i in range(len(standard_table)):
    fig = plt.figure(figsize=(15,3))
    continuum, region_indices = ig.local_continuum_fit(wavelen_stack[:,i],
                                                    raw_flux_stack[:,i],
                                                    poly_order = poly_deg,
                                                    line_center = line_center,
                                                    spec_res = spec_res,
                                                    regions = regions)
    # append indices to list of indices
    reg_idx_stack.append(region_indices)
    
    # append to list of the local continuum arrays
    continuum_stack.append(continuum)

    # normalize flux by dividing the flux by the continuum and create a list
    norm_flux = raw_flux_stack[region_indices[0][0]:region_indices[n][1],i]/continuum_stack[i]
    norm_flux_stack.append(norm_flux)


    plt.plot(wavelen_stack[region_indices[0][0]:region_indices[n][1],i], raw_flux_stack[region_indices[0][0]:region_indices[n][1],i])
    plt.plot(wavelen_stack[region_indices[0][0]:region_indices[n][1],i], continuum_stack[i])

    for j in range(len(region_indices)):
            plt.axvspan(wavelen_stack[region_indices[n-j][0],i],wavelen_stack[region_indices[n-j][1],i],color='red',alpha=.2)

    plt.title(f"{line_name} {standard_table['Name'][i]} {standard_table['Spectral_Type'][i]}")

    plt.show()
    plt.show()

# Initialize storage lists
params_stack = []
params_error_stack = []
best_model_stack = []
result_stack = []
dely_stack = []

# Define initial parameters for Gaussian fitting
init_params = (-0.5, line_center, spec_res,
               -0.5, lines_table[al1_mask]['Observed'][0], spec_res)

for i in range(len(standard_table)):    
    # Define initial parameters for Gaussian fitting
    params = Parameters()

    params.add('amp1', value = 1e-5, max = 0)
    params.add('c1', value = line_center)
    params.add('std1', value = spec_res)

    params.add('amp2', value = 1e-5, max = 0)
    params.add('c2', value = lines_table[al1_mask]['Observed'][0])
    params.add('std2', value = spec_res)

    # params['std1'].expr = 'std2' # constrain widths to be the same

    # Define Model
    model = Model(ig.two_gaussian, nan_policy='omit')

    result = model.fit(norm_flux_stack[i], params,
                       x = wavelen_stack[reg_idx_stack[i][0][0]:reg_idx_stack[i][n][1], i],
                       weights = 1/raw_flux_err_stack[reg_idx_stack[i][0][0]:reg_idx_stack[i][n][1], i])
    result_stack.append(result) # lmfit ModelResult objects
    dely = result.eval_uncertainty(sigma=5) # lmfit
    dely_stack.append(dely)

    # Append the best_model, error and parameters stacks
    # could probably make a loop to the tune of for value in result etc...
    params_stack.extend([result.params['amp1'].value, result.params['c1'].value, result.params['std1'].value,
                         result.params['amp2'].value, result.params['c2'].value, result.params['std2'].value])
    
    params_error_stack.extend([result.params['amp1'].stderr, result.params['c1'].stderr, result.params['std1'].stderr,
                               result.params['amp2'].stderr, result.params['c2'].stderr, result.params['std2'].stderr])

    best_model_stack.append(result.best_fit)

params_arr = np.array(params_stack)

amps = params_arr[0::3]
amps_err = params_error_stack[0::3]

centers = params_arr[1::3]
center_err = params_error_stack[1::3]

sigmas = params_arr[2::3]
sigmas_err = params_error_stack[2::3]

# Gaussian model fits for each source
mg1_2_1065_fits = []
al1_2_1098_fits = []
num_gauss = 2

# Getting the individual Gaussians into their own lists
for i in range(len(standard_table)):
    mg1_2_1065_fits.append(ig.gaussian(wavelen_stack[reg_idx_stack[i][0][0]:reg_idx_stack[i][n][1],i],amps[0::num_gauss][i],centers[0::num_gauss][i],sigmas[0::num_gauss][i]))
    al1_2_1098_fits.append(ig.gaussian(wavelen_stack[reg_idx_stack[i][0][0]:reg_idx_stack[i][n][1],i],amps[1::num_gauss][i],centers[1::num_gauss][i],sigmas[1::num_gauss][i]))

ew1_stack = []
ew2_stack = []
# ew3_stack = []
# ew4_stack = []
for i in range(len(standard_table)):
    ew1 = np.trapz(1-mg1_2_1065_fits[i], wavelen_stack[reg_idx_stack[i][0][0]:reg_idx_stack[i][n][1],i])
    ew2 = np.trapz(1-al1_2_1098_fits[i], wavelen_stack[reg_idx_stack[i][0][0]:reg_idx_stack[i][n][1],i])

    ew1_stack.append(ew1)
    ew2_stack.append(ew2)

standard_table['ew_mg1_2_1065'] = ew1_stack
standard_table['ew_al1_2_1098'] = ew2_stack

# Plotting each component Gaussian and the best model fit
for i in range(len(standard_table)):
    # Create subplots with adjusted size ratios
    fig, (ax1, ax2) = plt.subplots(2, 1, figsize=(15, 6), gridspec_kw={'height_ratios': [4, 1]}, sharex=True)

    ax1.plot(wavelen_stack[reg_idx_stack[i][0][0]:reg_idx_stack[i][n][1],i], best_model_stack[i], c='red', label='Gaussian Fit')
    ax1.plot(wavelen_stack[reg_idx_stack[i][0][0]:reg_idx_stack[i][n][1],i], norm_flux_stack[i], c='black', alpha=0.5)
    ax1.fill_between(wavelen_stack[reg_idx_stack[i][0][0]:reg_idx_stack[i][n][1],i],
                     best_model_stack[i]-dely_stack[i], best_model_stack[i]+dely_stack[i],
                     alpha=0.2, color='red')

    ax1.plot(wavelen_stack[reg_idx_stack[i][0][0]:reg_idx_stack[i][n][1],i], mg1_2_1065_fits[i], ls=':',
     label=r'Mg I 2.1065 $\mu$m' '\n' rf'A={amps[0::num_gauss][i]:.3e}' '\n' rf'$\mu$={centers[0::num_gauss][i]:.5f}' '\n' rf'$\sigma$={sigmas[0::num_gauss][i]:.3e}')
    ax1.plot(wavelen_stack[reg_idx_stack[i][0][0]:reg_idx_stack[i][n][1],i], al1_2_1098_fits[i], ls=':',
     label=r'Al I 2.1098 $\mu$m' '\n' rf'A={amps[1::num_gauss][i]:.3e}' '\n' rf'$\mu$={centers[1::num_gauss][i]:.5f}' '\n' rf'$\sigma$={sigmas[1::num_gauss][i]:.3e}')

    ax1.set_ylabel('Flux + Constant')
    ax1.legend(bbox_to_anchor=(1,1))

    # Calculate residuals
    residuals = (norm_flux_stack[i] - best_model_stack[i])
    sqsum_res = np.sum(residuals**2)
    sqsum = np.sum((norm_flux_stack[i]-np.mean(norm_flux_stack[i])**2))
    R2 = 1-(sqsum_res/sqsum)

    # Plot residuals
    ax2.scatter(wavelen_stack[reg_idx_stack[i][0][0]:reg_idx_stack[i][n][1], i], residuals, c='red', label='Residuals', s=3)
    ax2.axhline(y=0, color='grey', linestyle='--')
    # ax2.set_ylim(-0.1,0.1)
    ax2.set_xlabel(r'Wavelength [$\mu$m]')
    ax2.set_ylabel('res')
    ax2.set_title(rf'$R^2$={R2:.3f}', fontsize=12)

    plt.suptitle(rf"Mg/Al Interval {standard_table['Name'][i]} {standard_table['Spectral_Type'][i]}")
    plt.show()

fig = plt.figure(figsize=(15,5))

plt.plot(standard_table['Spectral_Type'], standard_table['ew_mg1_2_1065'],ls='--',marker='p', label=f"Mg I 2.1065 $\mu$m")
plt.plot(standard_table['Spectral_Type'], standard_table['ew_al1_2_1098'],ls='--',marker='p', label=f"Al I 2.1098 $\mu$m")
# plt.plot(standard_table['Spectral_Type'], standard_table['ew_fe1_2_2838'],ls='--',marker='p', label=f"Fe I 2.2838$\mu$m")

# plt.scatter(standard_table['Spectral_Type'], area_stack)
# plt.ylim(np.nanmedian(ew_stack)+-0.001,np.nanmedian(ew_stack)+0.001)

# plt.ylim(-0.0001,0.0002)
# plt.yscale('log')

plt.xlabel('Spectral Type')
plt.ylabel(r'Equivalent Width [$\mu$m]')
# plt.title(f"{line_name} {line_center} $\mu$m")

plt.legend(bbox_to_anchor=(1,1))
plt.show()

Ti Region (+ Fe)

# Define the region for fitting
line_name = lines_table[fe1_mask][0]['Spectrum'] # Species
line_center = lines_table[fe1_mask][85]['Observed'] # Wavelength

# from igrins_mod import local_continuum_fit
continuum_stack = []
norm_flux_stack = []
reg_idx_stack = []

# regions [left point,width]
regions = [(-600,30),(-350,50),(-150,60),(130,10),(400,100)]

# number of regions I use and subtract 1 since index starts at 0
n = len(regions)-1
poly_deg = 3
for i in range(len(standard_table)):
    fig = plt.figure(figsize=(15,3))
    continuum, region_indices = ig.local_continuum_fit(wavelen_stack[:,i],
                                                    raw_flux_stack[:,i],
                                                    poly_order = poly_deg,
                                                    line_center = line_center,
                                                    spec_res = spec_res,
                                                    regions = regions)
    # append indices to list of indices
    reg_idx_stack.append(region_indices)
    
    # append to list of the local continuum arrays
    continuum_stack.append(continuum)

    # normalize flux by dividing the flux by the continuum and create a list
    norm_flux = raw_flux_stack[region_indices[0][0]:region_indices[n][1],i]/continuum_stack[i]
    norm_flux_stack.append(norm_flux)


    plt.plot(wavelen_stack[region_indices[0][0]:region_indices[n][1],i], raw_flux_stack[region_indices[0][0]:region_indices[n][1],i])
    plt.plot(wavelen_stack[region_indices[0][0]:region_indices[n][1],i], continuum_stack[i])

    for j in range(len(region_indices)):
        plt.axvspan(wavelen_stack[region_indices[n-j][0],i],wavelen_stack[region_indices[n-j][1],i],color='red',alpha=.2)

    plt.title(f"{line_name} {standard_table['Name'][i]} {standard_table['Spectral_Type'][i]}")

    plt.show()

params_stack = []
params_error_stack = []

pcov_stack = []
best_model_stack = []

# Define initial parameters for Gaussian fitting
init_params = (-0.2, lines_table[ti1_mask]['Observed'][44], spec_res,
               -0.2, lines_table[ti1_mask]['Observed'][45], spec_res,
                -0.2, line_center, spec_res,
                -0.1, lines_table[fe1_mask]['Observed'][86], spec_res,
                -0.1, lines_table[ti1_mask]['Observed'][47], spec_res)

for i in range(len(standard_table)):    
    # Perform Gaussian fitting for the current source
    popt, pcov, param_err, best_model = ig.model_fit(ig.five_gaussian,
                                        wavelen_stack[reg_idx_stack[i][0][0]:reg_idx_stack[i][n][1],i],
                                        norm_flux_stack[i],
                                        raw_flux_err_stack[reg_idx_stack[i][0][0]:reg_idx_stack[i][n][1],i],
                                        init_params,
                                        maxfev=5000)
    # Append the best_model, error and parameters stacks
    params_error_stack.extend(param_err)

    params_stack.extend(popt)
    params_arr = np.array(params_stack)

    amps = params_arr[0::3]
    amps_err = params_error_stack[0::3]

    centers = params_arr[1::3]
    center_err = params_error_stack[1::3]

    sigmas = params_arr[2::3]
    sigmas_err = params_error_stack[2::3]

    pcov_stack.append(pcov)
    best_model_stack.append(best_model)

ti1_2_2217_fit = []
ti1_2_2239_fit = []
fe1_2_2263_fit = [] # maybe fix wavelength offsets
fe1_2_2266_fit = []
ti1_2_2280_fit = []
num_gauss = 5

for i in range(len(standard_table)):
    ti1_2_2217_fit.append(ig.gaussian(wavelen_stack[reg_idx_stack[i][0][0]:reg_idx_stack[i][n][1],i],
        amps[0::num_gauss][i],centers[0::num_gauss][i],sigmas[0::num_gauss][i]))

    ti1_2_2239_fit.append(ig.gaussian(wavelen_stack[reg_idx_stack[i][0][0]:reg_idx_stack[i][n][1],i],
        amps[1::num_gauss][i],centers[1::num_gauss][i],sigmas[1::num_gauss][i]))

    fe1_2_2263_fit.append(ig.gaussian(wavelen_stack[reg_idx_stack[i][0][0]:reg_idx_stack[i][n][1],i],
        amps[2::num_gauss][i],centers[2::num_gauss][i],sigmas[2::num_gauss][i]))

    fe1_2_2266_fit.append(ig.gaussian(wavelen_stack[reg_idx_stack[i][0][0]:reg_idx_stack[i][n][1],i],
        amps[3::num_gauss][i],centers[3::num_gauss][i],sigmas[3::num_gauss][i]))

    ti1_2_2280_fit.append(ig.gaussian(wavelen_stack[reg_idx_stack[i][0][0]:reg_idx_stack[i][n][1],i],
        amps[4::num_gauss][i],centers[4::num_gauss][i],sigmas[4::num_gauss][i]))
    

# Integate 1-best_model to get the areas (equivalent widths) of each component Gaussian
for i in range(len(standard_table)):
    # Create subplots with adjusted size ratios
    fig, (ax1, ax2) = plt.subplots(2, 1, figsize=(15, 6), gridspec_kw={'height_ratios': [4, 1]}, sharex=True)

    ax1.plot(wavelen_stack[reg_idx_stack[i][0][0]:reg_idx_stack[i][n][1],i], norm_flux_stack[i], alpha=0.2, c='black', label='Data')
    ax1.plot(wavelen_stack[reg_idx_stack[i][0][0]:reg_idx_stack[i][n][1],i], best_model_stack[i], c='red', lw=2, label='Best Model')
    
    ax1.plot(wavelen_stack[reg_idx_stack[i][0][0]:reg_idx_stack[i][n][1],i],ti1_2_2217_fit[i], ls='--',
        label = 'Ti I 2.2217 Fit' '\n' rf'A={amps[0::num_gauss][i]:.3e}' '\n' rf'$\mu$={centers[0::num_gauss][i]:.5f}' '\n' rf'$\sigma$={sigmas[0::num_gauss][i]:.3e}')

    ax1.plot(wavelen_stack[reg_idx_stack[i][0][0]:reg_idx_stack[i][n][1],i],ti1_2_2239_fit[i], ls='--',
        label = 'Ti I 2.2239 Fit' '\n' rf'A={amps[1::num_gauss][i]:.3e}' '\n' rf'$\mu$={centers[1::num_gauss][i]:.5f}' '\n' rf'$\sigma$={sigmas[1::num_gauss][i]:.3e}')

    ax1.plot(wavelen_stack[reg_idx_stack[i][0][0]:reg_idx_stack[i][n][1],i],fe1_2_2263_fit[i], ls='--',
        label = 'Fe I 2.2263 Fit' '\n' rf'A={amps[2::num_gauss][i]:.3e}' '\n' rf'$\mu$={centers[2::num_gauss][i]:.5f}' '\n' rf'$\sigma$={sigmas[2::num_gauss][i]:.3e}')

    ax1.plot(wavelen_stack[reg_idx_stack[i][0][0]:reg_idx_stack[i][n][1],i],fe1_2_2266_fit[i], ls='--',
        label = 'Fe I 2.2266 Fit' '\n' rf'A={amps[3::num_gauss][i]:.3e}' '\n' rf'$\mu$={centers[3::num_gauss][i]:.5f}' '\n' rf'$\sigma$={sigmas[3::num_gauss][i]:.3e}')

    ax1.plot(wavelen_stack[reg_idx_stack[i][0][0]:reg_idx_stack[i][n][1],i],ti1_2_2280_fit[i], ls='--',
        label = 'Ti I 2.2280 Fit' '\n' rf'A={amps[4::num_gauss][i]:.3e}' '\n' rf'$\mu$={centers[4::num_gauss][i]:.5f}' '\n' rf'$\sigma$={sigmas[4::num_gauss][i]:.3e}')

 
    ax1.set_ylabel('Flux + Constant')
    ax1.legend(bbox_to_anchor=(1,1), fontsize=12)

    # Calculate residuals
    residuals = (norm_flux_stack[i] - best_model_stack[i])
    sqsum_res = np.sum(residuals**2)
    sqsum = np.sum((norm_flux_stack[i]-np.mean(norm_flux_stack[i])**2))
    R2 = 1-(sqsum_res/sqsum)

    # Plot residuals
    ax2.scatter(wavelen_stack[reg_idx_stack[i][0][0]:reg_idx_stack[i][n][1], i], residuals, c='red', label='Residuals', s=3)
    ax2.axhline(y=0, color='grey', linestyle='--')
    # ax2.set_ylim(-0.1,0.1)
    ax2.set_xlabel(r'Wavelength [$\mu$m]')
    ax2.set_ylabel('res')
    ax2.set_title(rf'$R^2$={R2:.3f}', fontsize=12)

    plt.suptitle(rf"Ti Interval {standard_table['Name'][i]} {standard_table['Spectral_Type'][i]}")

    plt.show()

ew1_stack = [] # empty list to load in equivalent widths
ew2_stack = []
ew3_stack = []

# Integate 1-best_model to get the areas (equivalent widths) of each component Gaussian
for i in range(len(standard_table)):
    ew1 = np.trapz(1-fe1_2_2263_fit[i],wavelen_stack[reg_idx_stack[i][0][0]:reg_idx_stack[i][n][1],i]) # integrate the gaussian
    ew1_stack.append(ew1)

    ew2 = np.trapz(1-fe1_2_2266_fit[i],wavelen_stack[reg_idx_stack[i][0][0]:reg_idx_stack[i][n][1],i])
    ew2_stack.append(ew2)

    ew3 = np.trapz(1-ti1_2_2280_fit[i],wavelen_stack[reg_idx_stack[i][0][0]:reg_idx_stack[i][n][1],i])
    ew3_stack.append(ew3)

standard_table.loc[:,'ew_fe1_2_2263'] = ew1_stack
standard_table.loc[:,'ew_fe1_2_2266'] = ew2_stack

standard_table.loc[:,'ew_ti1_2_2280'] = ew3_stack

fig = plt.figure(figsize=(15,5))

plt.plot(standard_table['Spectral_Type'], np.abs(standard_table['ew_fe1_2_2263']),ls='-',marker='o', label=f"Fe I 2.2263 $\mu$m")
plt.plot(standard_table['Spectral_Type'], np.abs(standard_table['ew_fe1_2_2266']),ls='-',marker='o', label=f"Fe I 2.2266 $\mu$m")
plt.plot(standard_table['Spectral_Type'], np.abs(standard_table['ew_ti1_2_2280']),ls='-',marker='o', label=f"Fe I 2.2280 $\mu$m")

plt.xlabel('Spectral Type')
plt.ylabel(r'Equivalent Width [$\mu$m]')
# plt.title(f"{line_name} {line_center} $\mu$m")


plt.legend(bbox_to_anchor=(1, 1))
plt.show()

Ca Region

Ca I 2.2614

Fe I 2.2626 = lines_table[fe1_mask][29]

Ca I 2.2631

Ca I 2.2657

# import matplotlib
# matplotlib.use('TkAgg')  # Set the backend to TkAgg for interactivity
# # Define the region for fitting
# line_name = lines_table[ca1_mask][0]['Spectrum'] # Species
# line_center = lines_table[ca1_mask][0]['Observed'] # Wavelength

# # from igrins_mod import local_continuum_fit
# continuum_stack = []
# norm_flux_stack = []
# reg_idx_stack = []

# # Helper function to select regions interactively
# def select_regions(wavelengths, flux, num_regions):
#     plt.plot(wavelengths, flux)
#     plt.title('Click to select points for region exclusion')
#     plt.xlabel('Wavelength')
#     plt.ylabel('Flux')

#     # Select points interactively
#     selected_points = plt.ginput(n=num_regions*2, timeout=0) # n=num_regions*2 because we need start and end points for each region
#     plt.close()

#     regions = []
#     for i in range(0, len(selected_points), 2):
#         start = selected_points[i][0]
#         end = selected_points[i+1][0]
#         regions.append((start, end))
    
#     return regions

# # number of regions to select interactively
# num_regions = 4
# poly_deg = 3

# for i in range(len(standard_table)):
#     idx1 = np.searchsorted(wavelen_stack[:, i], line_center ) - 150
#     idx2 = np.searchsorted(wavelen_stack[:, i], lines_table[ca1_mask]['Observed'][2]) + 150

#     wavelengths = wavelen_stack[idx1:idx2, i]
#     flux = raw_flux_stack[idx1:idx2, i]
    
#     # Select regions interactively
#     regions = select_regions(wavelengths, flux, num_regions)
    
#     # Prepare the regions for the local_continuum_fit function
#     regions = [(int(np.searchsorted(wavelengths, start)), int(np.searchsorted(wavelengths, end))) for start, end in regions]
    
#     fig = plt.figure(figsize=(15, 3))
#     continuum, region_indices = ig.local_continuum_fit(wavelengths,
#                                                        flux,
#                                                        poly_order=poly_deg,
#                                                        line_center=line_center,
#                                                        spec_res=spec_res,
#                                                        regions=regions)
#     # Append indices to list of indices
#     reg_idx_stack.append(region_indices)
    
#     # Append to list of the local continuum arrays
#     continuum_stack.append(continuum)

#     # Normalize flux by dividing the flux by the continuum and create a list
#     norm_flux = flux[region_indices[0][0]:region_indices[-1][1]] / continuum_stack[i]
#     norm_flux_stack.append(norm_flux)

#     plt.plot(wavelengths[region_indices[0][0]:region_indices[-1][1]], flux[region_indices[0][0]:region_indices[-1][1]])
#     plt.plot(wavelengths[region_indices[0][0]:region_indices[-1][1]], continuum_stack[i])

#     for j in range(len(region_indices)):
#         plt.axvspan(wavelengths[region_indices[j][0]], wavelengths[region_indices[j][1]], color='red', alpha=.2)

#     plt.title(f"{line_name} {standard_table['Name'][i]} {standard_table['Spectral_Type'][i]}")

#     plt.show()

# Define the region for fitting
line_name = lines_table[ca1_mask][0]['Spectrum'] # Species
line_center = lines_table[ca1_mask][0]['Observed'] # Wavelength
# from igrins_mod import local_continuum_fit
continuum_stack = []
norm_flux_stack = []
reg_idx_stack = []

# regions [left point,width]
regions = [(-100,50), (50,30), (240,50),(550,10)]

# number of regions I use and subtract 1 since index starts at 0
n = len(regions)-1
poly_deg = 3
for i in range(len(standard_table)):
    fig = plt.figure(figsize=(15,3))
    continuum, region_indices = ig.local_continuum_fit(wavelen_stack[:,i],
                                                    raw_flux_stack[:,i],
                                                    poly_order = poly_deg,
                                                    line_center = line_center,
                                                    spec_res = spec_res,
                                                    regions = regions)
    # append indices to list of indices
    reg_idx_stack.append(region_indices)
    
    # append to list of the local continuum arrays
    continuum_stack.append(continuum)

    # normalize flux by dividing the flux by the continuum and create a list
    norm_flux = raw_flux_stack[region_indices[0][0]:region_indices[n][1],i]/continuum_stack[i]
    norm_flux_stack.append(norm_flux)


    plt.plot(wavelen_stack[region_indices[0][0]:region_indices[n][1],i], raw_flux_stack[region_indices[0][0]:region_indices[n][1],i])
    plt.plot(wavelen_stack[region_indices[0][0]:region_indices[n][1],i], continuum_stack[i])

    for j in range(len(region_indices)):
        plt.axvspan(wavelen_stack[region_indices[n-j][0],i],wavelen_stack[region_indices[n-j][1],i],color='red',alpha=.2)

    plt.title(f"{line_name} {standard_table['Name'][i]} {standard_table['Spectral_Type'][i]}")

    plt.show()

np.abs(lines_table[ca1_mask]['Observed'][0] - lines_table[ca1_mask]['Observed'][1])
# lines_table[ca1_mask]['Observed'][1]
# lines_table[ca1_mask]['Observed'][2]
# lines_table[fe1_mask]['Observed'][104]

0.0017000000000000348

# Initialize storage lists
params_stack = []
params_error_stack = []
best_model_stack = []
result_stack = []
dely_stack = []

# Define initial parameters for Gaussian fitting
# init_params = (-0.5, lines_table[ca1_mask]['Observed'][0], spec_res,
#                -0.5, lines_table[ca1_mask]['Observed'][1], spec_res,
#                -0.5, lines_table[ca1_mask]['Observed'][2], spec_res,
#                -0.5, lines_table[fe1_mask]['Observed'][104], spec_res)
from scipy.signal import savgol_filter



for i in range(len(standard_table)):    
    # Define initial parameters for Gaussian fitting
    params = Parameters()

    params.add('amp1', value = 1e-5, max = 0)
    params.add('c1', value = lines_table[ca1_mask]['Observed'][0])
    params.add('std1', value = spec_res)

    params.add('amp2', value = 1e-5, max = 0)
    params.add('c2', value = lines_table[ca1_mask]['Observed'][1])
    params.add('std2', value = spec_res)

    params.add('amp3', value = 1e-5, max = 0)
    params.add('c3', value = lines_table[ca1_mask]['Observed'][2])
    params.add('std3', value = spec_res)

    params.add('amp4', value = 1e-5, max = 0)
    params.add('c4', value = lines_table[fe1_mask]['Observed'][104])
    params.add('std4', value = spec_res)

    params['std1'].expr = 'std2' # constrain widths to be the same
    params['std1'].expr = 'std3'
    # params['std1'].expr = 'std2'

    # Define Model
    model = Model(ig.four_gaussian, nan_policy='omit')

    result = model.fit(savgol_filter(norm_flux_stack[i], window_length=11, polyorder=3), params,
                       x = wavelen_stack[reg_idx_stack[i][0][0]:reg_idx_stack[i][n][1], i],
                       weights = 1/raw_flux_err_stack[reg_idx_stack[i][0][0]:reg_idx_stack[i][n][1], i])

    result_stack.append(result) # lmfit ModelResult objects
    dely = result.eval_uncertainty(sigma=5) # lmfit
    dely_stack.append(dely)

    # Append the best_model, error and parameters stacks
    # could probably make a loop to the tune of for value in result etc...
    params_stack.extend([result.params['amp1'].value, result.params['c1'].value, result.params['std1'].value,
                         result.params['amp2'].value, result.params['c2'].value, result.params['std2'].value,
                         result.params['amp3'].value, result.params['c3'].value, result.params['std3'].value,
                         result.params['amp4'].value, result.params['c4'].value, result.params['std4'].value])
    
    params_error_stack.extend([result.params['amp1'].stderr, result.params['c1'].stderr, result.params['std1'].stderr,
                               result.params['amp2'].stderr, result.params['c2'].stderr, result.params['std2'].stderr,
                               result.params['amp3'].stderr, result.params['c3'].stderr, result.params['std3'].stderr,
                               result.params['amp4'].stderr, result.params['c4'].stderr, result.params['std4'].stderr])

    best_model_stack.append(result.best_fit)

params_arr = np.array(params_stack)

amps = params_arr[0::3]
amps_err = params_error_stack[0::3]

centers = params_arr[1::3]
center_err = params_error_stack[1::3]

sigmas = params_arr[2::3]
sigmas_err = params_error_stack[2::3]

# Gaussian model fits for each source
ca1_2_2614_fits = []
ca1_2_2631_fits = []
ca1_2_2657_fits = []
fe1_2_2626_fits = []
num_gauss = 4

# Getting the individual Gaussians into their own lists
for i in range(len(standard_table)):
    ca1_2_2614_fits.append(ig.gaussian(wavelen_stack[reg_idx_stack[i][0][0]:reg_idx_stack[i][n][1], i], amps[0::num_gauss][i], centers[0::num_gauss][i], sigmas[0::num_gauss][i]))
    ca1_2_2631_fits.append(ig.gaussian(wavelen_stack[reg_idx_stack[i][0][0]:reg_idx_stack[i][n][1], i], amps[1::num_gauss][i], centers[1::num_gauss][i], sigmas[1::num_gauss][i]))
    ca1_2_2657_fits.append(ig.gaussian(wavelen_stack[reg_idx_stack[i][0][0]:reg_idx_stack[i][n][1], i], amps[2::num_gauss][i], centers[2::num_gauss][i], sigmas[2::num_gauss][i]))
    fe1_2_2626_fits.append(ig.gaussian(wavelen_stack[reg_idx_stack[i][0][0]:reg_idx_stack[i][n][1], i], amps[3::num_gauss][i], centers[3::num_gauss][i], sigmas[3::num_gauss][i]))

# Plotting each component Gaussian and the best model fit
for i in range(len(standard_table)):
    # Create subplots with adjusted size ratios
    fig, (ax1, ax2) = plt.subplots(2, 1, figsize=(15, 6), gridspec_kw={'height_ratios': [4, 1]}, sharex=True)

    ax1.plot(wavelen_stack[reg_idx_stack[i][0][0]:reg_idx_stack[i][n][1],i], best_model_stack[i], c='red', label=r'Best Fit $\pm 5 \sigma$')
    ax1.plot(wavelen_stack[reg_idx_stack[i][0][0]:reg_idx_stack[i][n][1],i], savgol_filter(norm_flux_stack[i], window_length=11, polyorder=3), c='black', alpha=0.5)

    ax1.fill_between(wavelen_stack[reg_idx_stack[i][0][0]:reg_idx_stack[i][n][1],i],
                    best_model_stack[i]-dely_stack[i], best_model_stack[i]+dely_stack[i],
                    alpha=0.2, color = 'red' )

    ax1.plot(wavelen_stack[reg_idx_stack[i][0][0]:reg_idx_stack[i][n][1],i], ca1_2_2614_fits[i], ls='-.',
            label=r'Ca I 2.2614 $\mu$m' "\n" rf"A={amps[0::num_gauss][i]:.3e}" "\n" rf"$\mu$={centers[0::num_gauss][i]:.5f}" "\n" rf"$\sigma$= {sigmas[0::num_gauss][i]:.3e}")
    
    ax1.plot(wavelen_stack[reg_idx_stack[i][0][0]:reg_idx_stack[i][n][1],i], ca1_2_2631_fits[i], ls='-.',
        label=r'Ca I 2.2631 $\mu$m' "\n" rf"A={amps[0::num_gauss][i]:.3e}" "\n" rf"$\mu$={centers[0::num_gauss][i]:.5f}" "\n" rf"$\sigma$= {sigmas[0::num_gauss][i]:.3e}")
    ax1.plot(wavelen_stack[reg_idx_stack[i][0][0]:reg_idx_stack[i][n][1],i], ca1_2_2657_fits[i], ls='-.',
        label=r'Ca I 2.2657 $\mu$m' "\n" rf"A={amps[0::num_gauss][i]:.3e}" "\n" rf"$\mu$={centers[0::num_gauss][i]:.5f}" "\n" rf"$\sigma$= {sigmas[0::num_gauss][i]:.3e}")
    ax1.plot(wavelen_stack[reg_idx_stack[i][0][0]:reg_idx_stack[i][n][1],i], fe1_2_2626_fits[i], ls='-.', c='magenta',
        label=r'Fe I 2.2626 $\mu$m' "\n" rf"A={amps[0::num_gauss][i]:.3e}" "\n" rf"$\mu$={centers[0::num_gauss][i]:.5f}" "\n" rf"$\sigma$= {sigmas[0::num_gauss][i]:.3e}")
    
    # plt.axvline()
    plt.suptitle(rf"Ca Region {standard_table['Name'][i]} {standard_table['Spectral_Type'][i]}")

    ax1.set_ylabel('Flux + Constant')
    ax1.legend(bbox_to_anchor=(1,1))

    # Calculate residuals
    residuals = (norm_flux_stack[i] - best_model_stack[i])
    sqsum_res = np.sum(residuals**2)
    sqsum = np.sum((norm_flux_stack[i]-np.mean(norm_flux_stack[i])**2))
    R2 = 1-(sqsum_res/sqsum)

    # Plot residuals
    ax2.scatter(wavelen_stack[reg_idx_stack[i][0][0]:reg_idx_stack[i][n][1], i], residuals, c='red', label='Residuals', s=3)
    ax2.axhline(y=0, color='grey', linestyle='--')
    # ax2.set_ylim(-0.1,0.1)
    ax2.set_xlabel(r'Wavelength [$\mu$m]')
    ax2.set_ylabel('res')
    ax2.set_title(rf'$R^2$={R2:.3f}', fontsize=12)
    
    plt.show()

ew1_stack = []
ew2_stack = []
ew3_stack = []
ew4_stack = []

for i in range(len(standard_table)):
    ew1 = np.trapz(1-ca1_2_2614_fits[i],wavelen_stack[reg_idx_stack[i][0][0]:reg_idx_stack[i][n][1],i])
    ew2 = np.trapz(1-ca1_2_2631_fits[i],wavelen_stack[reg_idx_stack[i][0][0]:reg_idx_stack[i][n][1],i])
    ew3 = np.trapz(1-ca1_2_2657_fits[i],wavelen_stack[reg_idx_stack[i][0][0]:reg_idx_stack[i][n][1],i])
    ew4 = np.trapz(1-fe1_2_2626_fits[i],wavelen_stack[reg_idx_stack[i][0][0]:reg_idx_stack[i][n][1],i])

    ew1_stack.append(ew1)
    ew2_stack.append(ew2)
    ew3_stack.append(ew3)
    ew4_stack.append(ew4)

standard_table['ew_ca1_2_2614'] = ew1_stack
standard_table['ew_ca1_2_2631'] = ew2_stack
standard_table['ew_ca1_2_2657'] = ew3_stack
standard_table['ew_fe1_2_2626'] = ew4_stack

fig = plt.figure(figsize=(15,5))

# plt.plot(standard_table['Spectral_Type'], standard_table['ew_sc1_2_2058'],ls='--',marker='p', label=f"Sc I 2.2058 $\mu$m")
# plt.plot(standard_table['Spectral_Type'], standard_table['ew_na1_2_2062'],ls='--',marker='p', label=f"Na I 2.2062 $\mu$m")
# plt.plot(standard_table['Spectral_Type'], standard_table['ew_si1_2_2068'],ls='--',marker='p', label=f"Si I 2.2068 $\mu$m")
# plt.plot(standard_table['Spectral_Type'], standard_table['ew_na1_2_2089'],ls='--',marker='p', label=f"Na I 2.2089 $\mu$m")

# plt.plot(standard_table['Spectral_Type'], np.abs(standard_table['ew_fe1_2_2263']),ls='-.',marker='o', label=f"Fe I 2.2263 $\mu$m")
# plt.plot(standard_table['Spectral_Type'], standard_table['ew_fe1_2_2266'],ls='-.',marker='o', label=f"Fe I 2.2266 $\mu$m")

plt.plot(standard_table['Spectral_Type'], standard_table['ew_ca1_2_2614'],ls=':',marker='s', label=f"Ca I 2.2614 $\mu$m")
plt.plot(standard_table['Spectral_Type'], standard_table['ew_ca1_2_2631'],ls=':',marker='s', label=f"Ca I 2.2631 $\mu$m")
plt.plot(standard_table['Spectral_Type'], standard_table['ew_ca1_2_2657'],ls=':',marker='s', label=f"Ca I 2.2657 $\mu$m")
plt.plot(standard_table['Spectral_Type'], standard_table['ew_fe1_2_2626'],ls=':',marker='s', label=f"Fe I 2.2626 $\mu$m")

# plt.scatter(standard_table['Spectral_Type'], area_stack)
# plt.ylim(np.nanmedian(ew_stack)+-0.001,np.nanmedian(ew_stack)+0.001)

# plt.ylim(-0.0001,0.0005)
# plt.yscale('log')

plt.xlabel('Spectral Type')
plt.ylabel(r'Equivalent Width [$\mu$m]')
# plt.title(f"{line_name} {line_center} $\mu$m")

plt.legend(bbox_to_anchor=(1, 1))
plt.show()

Mg, Ca, Fe ~ 2.28 micron

Mg I 2.2813 = lines_table[mg1_mask][8]

Fe I 2.2838 = lines_table[fe1_mask][115] (possible 116)

Ca ? = 2.285?

# Define the region for fitting
line_name = lines_table[mg1_mask][0]['Spectrum'] # Species
line_center = lines_table[mg1_mask][8]['Observed'] # Wavelength

# from igrins_mod import local_continuum_fit
continuum_stack = []
norm_flux_stack = []
reg_idx_stack = []

# regions [left point,width]
regions = [(-150,10),(-100,10),(70,10),(300,10),(500,10)]

# number of regions I use and subtract 1 since index starts at 0
n = len(regions)-1
poly_deg = 3
for i in range(len(standard_table)):
    fig = plt.figure(figsize=(15,3))
    continuum, region_indices = ig.local_continuum_fit(wavelen_stack[:,i],
                                                    raw_flux_stack[:,i],
                                                    poly_order = poly_deg,
                                                    line_center = line_center,
                                                    spec_res = spec_res,
                                                    regions = regions)
    # append indices to list of indices
    reg_idx_stack.append(region_indices)
    
    # append to list of the local continuum arrays
    continuum_stack.append(continuum)

    # normalize flux by dividing the flux by the continuum and create a list
    norm_flux = raw_flux_stack[region_indices[0][0]:region_indices[n][1],i]/continuum_stack[i]
    norm_flux_stack.append(norm_flux)


    plt.plot(wavelen_stack[region_indices[0][0]:region_indices[n][1],i], raw_flux_stack[region_indices[0][0]:region_indices[n][1],i])
    plt.plot(wavelen_stack[region_indices[0][0]:region_indices[n][1],i], continuum_stack[i])

    for j in range(len(region_indices)):
        plt.axvspan(wavelen_stack[region_indices[n-j][0],i],wavelen_stack[region_indices[n-j][1],i],color='red',alpha=.2)

    plt.title(f"{line_name} {standard_table['Name'][i]} {standard_table['Spectral_Type'][i]}")

    plt.show()

# min_wave = line_center - (100 * spec_res)
# max_wave = line_center + (100 * spec_res)

# Optimal parameters that will be appended from scipy curve fit below
amplitude1_stack = []
center1_stack = []
sigma1_stack = []

amplitude2_stack = []
centeR2_stack = []
sigma2_stack = []

amplitude3_stack = []
center3_stack = []
sigma3_stack = []

pcov_stack = [] # Covariant matrices
best_model_stack = [] # Best model fits

flux_constant = np.linspace(0,-1,len(standard_list))

# Define initial parameters for Gaussian fitting
init_params = ( -0.5, line_center, spec_res,
                -0.5, lines_table[fe1_mask]['Observed'][114], spec_res,
                -0.5, 2.2825, spec_res)

for i in range(len(standard_table)):    
    # Perform Gaussian fitting for the current source
    popt, pcov, param_err, best_model = ig.model_fit(ig.three_gaussian,
                                        wavelen_stack[reg_idx_stack[i][0][0]:reg_idx_stack[i][n][1],i],
                                        norm_flux_stack[i], raw_flux_err_stack[reg_idx_stack[i][0][0]:reg_idx_stack[i][n][1],i],
                                        init_params,
                                        maxfev=5000)
    # Append the Gaussian parameters, covariant matrices, and model fits as we loop through
    amplitude1_stack.append(popt[0])
    center1_stack.append(popt[1])
    sigma1_stack.append(popt[2])

    amplitude2_stack.append(popt[3])
    centeR2_stack.append(popt[4])
    sigma2_stack.append(popt[5])

    amplitude3_stack.append(popt[6])
    center3_stack.append(popt[7])
    sigma3_stack.append(popt[8])

    # amplitude4_stack.append(popt[7])
    # center4_stack.append(popt[8])

    pcov_stack.append(pcov)
    best_model_stack.append(best_model)

# Gaussian model fits for each source
mg1_2_2813_fits = []
fe1_2_2818_fits = []
ca1_2_2825_fits = []

# Getting the individual Gaussians into their own lists
for i in range(len(standard_table)):
    mg1_2_2813_fits.append(ig.gaussian(wavelen_stack[reg_idx_stack[i][0][0]:reg_idx_stack[i][n][1],i],amplitude1_stack[i],center1_stack[i],sigma1_stack[i]))
    fe1_2_2818_fits.append(ig.gaussian(wavelen_stack[reg_idx_stack[i][0][0]:reg_idx_stack[i][n][1],i],amplitude2_stack[i],centeR2_stack[i],sigma2_stack[i]))
    ca1_2_2825_fits.append(ig.gaussian(wavelen_stack[reg_idx_stack[i][0][0]:reg_idx_stack[i][n][1],i],amplitude3_stack[i],center3_stack[i],sigma3_stack[i]))

ew1_stack = []
ew2_stack = []
ew3_stack = []
# ew4_stack = []
for i in range(len(standard_table)):
    ew1 = np.trapz(1-mg1_2_2813_fits[i],wavelen_stack[reg_idx_stack[i][0][0]:reg_idx_stack[i][n][1],i])
    ew2 = np.trapz(1-fe1_2_2818_fits[i],wavelen_stack[reg_idx_stack[i][0][0]:reg_idx_stack[i][n][1],i])
    ew3 = np.trapz(1-ca1_2_2825_fits[i],wavelen_stack[reg_idx_stack[i][0][0]:reg_idx_stack[i][n][1],i])

    ew1_stack.append(ew1)
    ew2_stack.append(ew2)
    ew3_stack.append(ew3)


standard_table['ew_mg1_2_2813'] = ew1_stack
standard_table['ew_fe1_2_2818'] = ew2_stack
standard_table['ew_ca1_2_2825'] = ew3_stack

# Plotting each component Gaussian and the best model fit
for i in range(len(standard_table)):
    # Create subplots with adjusted size ratios
    fig, (ax1, ax2) = plt.subplots(2, 1, figsize=(15, 6), gridspec_kw={'height_ratios': [4, 1]}, sharex=True)

    ax1.plot(wavelen_stack[reg_idx_stack[i][0][0]:reg_idx_stack[i][n][1],i], best_model_stack[i], c='red',label='Gaussian Fit')
    ax1.plot(wavelen_stack[reg_idx_stack[i][0][0]:reg_idx_stack[i][n][1],i], norm_flux_stack[i], c='black', alpha=0.5)

    ax1.plot(wavelen_stack[reg_idx_stack[i][0][0]:reg_idx_stack[i][n][1],i], mg1_2_2813_fits[i],ls=':',label=r'Mg I 2.2813 $\mu$m')
    ax1.plot(wavelen_stack[reg_idx_stack[i][0][0]:reg_idx_stack[i][n][1],i], fe1_2_2818_fits[i],ls=':',label=r'Fe I 2.2818 $\mu$m')
    ax1.plot(wavelen_stack[reg_idx_stack[i][0][0]:reg_idx_stack[i][n][1],i], ca1_2_2825_fits[i],ls=':',label=r'Ca ? 2.2825 $\mu$m')
    
    # plt.axvline()
    ax1.suptitle(rf"Mg, Ca?, Fe {standard_table['Name'][i]} {standard_table['Spectral_Type'][i]}")

    ax1.set_ylabel('Flux + Constant')
    ax1.set_xlabel(r'Wavelength [$\mu$m]')
    ax1.legend(bbox_to_anchor=(1,1))
    
    plt.show()

c:\Users\Savio\anaconda3\envs\muler_dev\lib\site-packages\scipy\optimize\_minpack_py.py:1010: OptimizeWarning: Covariance of the parameters could not be estimated
  warnings.warn('Covariance of the parameters could not be estimated',

---------------------------------------------------------------------------
AttributeError                            Traceback (most recent call last)
Cell In[40], line 95
     92 ax1.plot(wavelen_stack[reg_idx_stack[i][0][0]:reg_idx_stack[i][n][1],i], ca1_2_2825_fits[i],ls=':',label=r'Ca ? 2.2825 $\mu$m')
     94 # plt.axvline()
---> 95 ax1.suptitle(rf"Mg, Ca?, Fe {standard_table['Name'][i]} {standard_table['Spectral_Type'][i]}")
     97 ax1.set_ylabel('Flux + Constant')
     98 ax1.set_xlabel(r'Wavelength [$\mu$m]')

AttributeError: 'Axes' object has no attribute 'suptitle'

fig = plt.figure(figsize=(15,5))

plt.plot(standard_table['Spectral_Type'], standard_table['ew_mg1_2_2813'],ls='--',marker='p', label=f"Mg I 2.2813 $\mu$m")
plt.plot(standard_table['Spectral_Type'], standard_table['ew_fe1_2_2818'],ls='--',marker='p', label=f"Fe I 2.2818 $\mu$m")
plt.plot(standard_table['Spectral_Type'], standard_table['ew_ca1_2_2825'],ls='--',marker='p', label=f"Ca I 2.2825 $\mu$m")

# plt.scatter(standard_table['Spectral_Type'], area_stack)
# plt.ylim(np.nanmedian(ew_stack)+-0.001,np.nanmedian(ew_stack)+0.001)

plt.ylim(-0.0001,0.0002)
# plt.yscale('log')

plt.xlabel('Spectral Type')
plt.ylabel(r'Equivalent Width [$\mu$m]')
# plt.title(f"{line_name} {line_center} $\mu$m")

plt.legend(bbox_to_anchor=(1,1))
plt.show()

# standard_table.to_csv('standard_table_v2.txt')